load(fullfile(matlabroot,'examples','signal','spots_num.mat'))

[pks,locs] = findpeaks(avSpots);
figure(1)
subplot(211);
plot(year,avSpots,year(locs),pks,'or')
xlabel Year, ylabel Number, axis tight

[pks,locs] = findpeaks(avSpots,'MinPeakDistance',6);
figure(1)
subplot(212);
plot(year,avSpots,year(locs),pks,'or')
xlabel Year, ylabel Number, title Sunspots, axis tight
legend('Data','peaks','location','NorthWest')

cycles = diff(locs); meanCycle = mean(cycles)
Fs = 1; Nf = 512;

df = Fs/Nf; f = 0:df:Fs/2-df;

trSpots = fftshift(fft(avSpots-mean(avSpots),Nf));

dBspots = 20*log10(abs(trSpots(Nf/2+1:Nf)));

yaxis = [20 85];
figure(2)
plot(f,dBspots,1./[meanCycle meanCycle],yaxis)
xlabel 'Frequency (year^{-1})', ylabel '| FFT | (dB)'
axis([0 1/2 yaxis])
text(1/meanCycle+.02,25,['<== 1/' num2str(meanCycle)])
[pk,MaxFreq] = findpeaks(dBspots,'NPeaks',1,'SortStr','descend');

Period = 1/f(MaxFreq)

hold on, plot(f(MaxFreq),pk,'or'), hold off
legend('Fourier transform','1/meanCycle','1/Period')